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Analytical and numerical calculations are presented for the mechanical response of fiber networks 
in a state of axisymmetric prestress, in the limit where geometric non-linearities such as fiber 
rotation are negligible. This allows us to focus on the anisotropy deriving purely from the non- 
linear force-extension curves of individual fibers. The number of independent elastic coefficients for 
isotropic, axisymmetric and fully anisotropic networks are enumerated, before deriving expressions 
for the response to a locally applied force that can be tested against e.g. microrheology experiments. 
Localised forces can generate anisotropy away from the point of application, so numerical integration 
of non-linear continuum equations is employed to determine the stress field, and induced mechanical 
anisotropy, at points located directly behind and in front of a force monopole. Results are presented 
for the wormlike chain model in normalised forms, allowing them to be easily mapped to a range of 
systems. Finally, the relevance of these findings to naturally occurring systems and directions for 
future investigation are discussed. 

PACS numbers: 87.16.Ka, 46.15.Cc, 46.25.Cc 



I. INTRODUCTION 

Many materials of industrial and biological importance 
are fibrous in nature, including paper pQ, carbon nan- 
otube assemblies [2J, and a range of protein fiber net- 
works such as the eukaryotic cellular cytoskeleton [3HS], 
some pathological and functional amyloids [BHH] j and self- 
assembled peptide networks [10H13] with applications in- 
cluding tissue engineering |14H16j . There now exists a 
substantial literature relating the macroscopic viscoelas- 
tic properties of such networks to their underlying mi- 
croscopic architecture [17l - f20] that has been verified for 
actin networks in particular [3TJ [55] , but also other pro- 
tein fibers such as vimentin ]23j . This theoretical frame- 
work is however not exhaustive, and noticeably the issue 
of anisotropy has received little attention to date, despite 
visualization of the actin cortex frequently demonstrat- 
ing a preferred orientation [H [SJ I24j . This morphologi- 
cal anisotropy has been been investigated in the context 
of coarse graining, highlighting relevant perturbations to 
the isotropic case |25j . and coupling to the environment 
leading to the alignment of stress fibers [2"u] . 

However, there is another form of mechanical 
anisotropy that arises, not from the geometric mi- 
crostructure of the network, but rather from the non- 
linear response of individual fibers. Consider applying 
an anisotropic prestress to an isotropic fiber network. If 
the magnitude of the stress is sufficient to place some 
fraction of fibers into the non-linear regime of their force- 
extension curves, the stiffness of individual fibers with 
respect to perturbations about this prestressed state will 
depend on their orientation, as schematically represented 
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in Fig. [T] The material response will therefore become 
anisotropic. Typically this prestress will also induce fiber 
rotation, but there is an important class of fiber networks 
for which this induced geometrical anisotropy can be ar- 
gued to be small. Many protein fibers have been found 
to be well described by the wormlike chain model, in 
which changes to the equilibrium end-to-end separation 
of fiber nodes induces an entropic restoring force |17j . For 
physiologically relevant parameters, such networks have 
been shown to strongly strain stiffen for modest strains of 
around 5-20% The geometrical anisotropy will thus 
remain small, even though the mechanical anisotropy is 
significant. This corresponds to geometrically linear elas- 
ticity with a non-linear constitutive equation, also known 
as hypoelastic elasticity (as opposed to hyperelastic elas- 
ticity where in addition the strains arc finite) 27 . 

The purpose of this article is to describe analytical and 
numerical calculations that quantitatively predict the 
mechanical response of fiber networks that have become 
anisotropic due to an axisymmetric prestress or prestrain. 
Changes to network geometry are entirely neglected, al- 
lowing the consequences of this form of anisotropy to be 
highlighted, while still generating results of relevance to 
many fiber networks. The primary assumptions are hy- 
poelasticity as described above, and also affinity, i.e. the 
strain field on fiber length scales is just a scaled-down 
version of the corresponding macroscopic strain. This 
strong assumption allows us to easily bridge the discrete 
and continuum representations. We also assume quasi- 
staticity, i.e. all calculations correspond to the elastic 
plateau regime of the network in question [T7] [TS] . 

Two key results are presented. Firstly, response func- 
tions relating the displacement caused by a locally ap- 
plied force are derived and reduced to a form that can be 
quickly integrated numerically, given network parameters 
and a prestress or prestrain. This can be employed in e.g. 
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active microrheology experiments, where a probe parti- 
cle is perturbed using an optical or magnetic trap [28] . 
allowing unknown parameters to be extracted via curve 
fitting to data. Secondly, the spatial stress field due to 
a force monopole is derived using a numerical procedure, 
and this is combined with the first result to derive predic- 
tions for the mechanical anisotropy induced by the force 
at a location directly in front or behind. This predictive 
procedure can again be employed to fit microrheology 
data, and is also relevant to networks with naturally- 
occurring force generators present, such as molecular mo- 
tors in some fiber protein networks. 



II. THEORY AND NUMERICAL METHOD 

In this section, analytical results for the mechanics of 
axisymmetric fiber networks are derived, including the 
local response to a point force. The numerical procedure 
for deriving the non-local stress field due to a point force 
is also explained. The notation used throughout is as 
follows. The line density (length of fibers per unit vol- 
ume) is denoted p, the persistence length of the fibers 
by £ p , and the distance between crosslinks by £ c . Where 
relevant, fiber orientation is denoted by the unit vector n. 



This paper is arranged as follows. In Sec. [TT] the key 
analytical results and numerical procedures are described 
in detail. The basic equations are presented in Sec. |II A| 
along with an enumeration of the reduction of indepen- 
dent elastic coefficients due to the specific microscopic 
picture assumed. Equations for the local response due 
to a point force are derived in Sec. |II B| for axisymmet- 
ric prestrain or prestress. The numerical procedure for 
deriving the spatial stress field in response to a localised 
force monopole is explained in Sec. |IIC[ Applications are 
presented in Sec. |III| where the wormlike chain model 
has been used throughout, allowing results to be pre- 
sented in normalised forms common to all fibers that 
obey this model. The elastic coefficients as a function 
of anisotropy and magnitude of either prestrain or pre- 
stress are presented in Sec. Ill A and the corresponding 



local response functions described in Sec. MB The me- 
chanical anisotropy induced by a localised force is given 
in Sec. |III C| Finally, the relevance of our results to real 
networks and future directions are discussed in Sec. Hvl 
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FIG. 1: Schematic of prestress-induced anisotropy. (a) The 
axes of principal stresses define Cartesian coordinates with the 
2-axis aligned with the maximum stress, (b) Fibers aligned 
parallel to this axis are placed in a tension r". (c) Those in 
the transverse direction have a tension r x , where t ± < t" . If 
the force-extension curve is non-linear, the network will thus 
become mechanically anisotropic, while potentially remaining 
(to good approximation) geometrically isotropic. 



A. Anisotropic linear network elasticity 

The linear constitutive equations for fiber networks un- 
dergoing affine deformation are now well established |17[ 
|2"9] . The assumptions here are that the unstrained net- 
work is isotropic, i.e. there is no net fiber orientation n, 
and that the strain uy remains small throughout. The 
stress Oij is then given by 



p{h i n ] T(£ c n k niu k i)) . 



(1) 



Here 5£ c = £ c n k niu k i is the extension of a fiber segment 
of length £ c and orientation n in the small-strain limit, 
t(5£ c ) is the corresponding tension, and (• • • )„ denotes 
averaging over all orientations n. ([I]) can be derived by 
e.g. considering the total vector force F across an arbi- 
trary plane with unit normal s, and using the definition 
of the stress tensor CijSj — Fi. 

The differential response to small changes in the strain 
uij + Suij follows from expanding ([I]) to linear order in 
Suij, giving a change in stress 

Saij = £ c p(T / (6£ c )h i n j n k n [ ) ii Su k i = C ijM 5u k i , (2) 

which defines the fourth-rank elastic stiffness ten- 
sor Cijki- Contributions due to fiber rotation are as- 
sumed to be sub-dominant and are ignored. In j2j, 
t'(5£ c ) denotes the gradient of the force-extension curve 
at the prestrained fiber extension S£ c , which depends on 
its orientation. If there is no prestrain, i.e. Uij = 0, it is 
customary to drop the <5's from ([2]), giving a more familiar 
expression relating Oij to uy. However, ^ still holds for 
perturbations about pre-stressed states uy ^ if the as- 
sumption of small strains remains valid, even though the 
force-extension relation t(S£ c ) may become non-linear. 
This corresponds to the hypoelastic case described in the 
introduction. 

By its definition in pb, the stiffness tensor CV,m is 
symmetric in all perturbations of its indices. This au- 
tomatically includes all of the symmetries required of 
any elastic stiffness tensor, i.e. Cij k i — Cji k i — C k iji — 
Ckiij OM [SB- It has 15 independent coefficients, which 
we can choose to be (in Cartesian coordinates) C xxxx , 

yyyy > ^zzzz: ^xxxy: ^xxxz, ^xxmi, ^xxzz, ^ XX nz, 



Cxyzzi ^xyyyi ^x 



cxxzi ^xxyyj ^xxzzt 



Cyyyz 5 ^yyZZ &nd Cy ZZZ . 



xxyz-j ^xyyzj 

This is a 



reduction from the 21 independent coefficients admitted 
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by an anisotropic elastic body in general [301 E] j which 
is due to the specific microscopic picture employed in de- 
riving ([2]). A similar reduction arises for axisymmetric 
systems as discussed in Sec. |II B| In fully isotropic sys- 
tems, where the Poisson ratio is fixed at 1/4 [321 [S3], there 
is only one independent coefficient (i.e. the shear modu- 
lus) rather than the usual two. It should be stressed this 
assumes affinity, without which ^ is not generally valid 
and the number of independent moduli may increase. 

B. Local mechanical response with axisymmetry 

We now turn to consider axisymmetric (also known as 
transversely isotropic |31| ) systems that are rotationally 
symmetric about a fixed axis. Without loss of generality, 
Cartesian coordinates (x, y, z) are used with the z-axis 
taken to be the axis of symmetry. The elastic stiffness 
tensor C^-jy now has fewer independent coefficients than 
the generally anisotropic case: Cijki is invariant under 
changes x -H> y, Cijki = if there are an odd number of 
x's or y's in the indices, and C xxxx — 3C xxyy as seen from 
performing the part of the integral ^ around the z-axis. 
Taking these additional constraints into account, there 



remain 3 independent elastic coefficients, here taken to 
be C xxxx = C , C xxzz = C 2 and C zzzz = C 4 . This is a re- 
duction from the 5 independent moduli for elastic bodies 
with the same symmetry |31j , again reflecting the specific 
microscopic picture assumed. A similar reduction (of 5 
to 3) occurs for nematic liquid crystals but for different 
reasons [31] , 

To determine the local mechanical response, a virtual 
point force f vlrt is applied at the origin, giving the force 
balance equations [30] . 

djSvy + /.r irt <5(x) = , (3) 

where the delta function localises the force at the origin. 
The djSdij in (|3| can be written in terms of the first 
derivative of the strain fluctuations 8v,ij using ^ (as- 
suming constant C^ki), and hence the second-derivatives 
of the displacement fluctuations Sui using standard for- 
mulae [3D]. This is then Fourier transformed as per 
<5{ii(q) = J e~ lq x <5zii(x)dx to give the matrix equation 

f virt = MSu (4) 

where using the notation q = {q x ,q y ,q z ), 
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Coq x + |Co9y 4 
^Coq x q y 
2C 2 q x q z 



C 2 q 2 



\Coq x q y 



\C q 2 x + C ql + C 2 q 2 z 
2C 2 q y q z 



2C 2 q x q z 
2C 2 q y q z 
C 2 (q 2 x + ql) + C iQ 2 z 



In order to obtain the spatial response field due to f vlrt ; 
M is inverted and inserted into Q, which would properly 
then be inverse Fourier transformed to give (5u(x); how- 
ever this inverse transform is not straightforward to per- 
form. Instead we restrict attention to the displacement 

u probc = ( M PiobO ; u probO ; u probo) of & sphere Q f radius a &t 

the origin to which the force f virt is applied, which allows 
an approximate solution to be readily attained [35j . In 



this procedure, the q-modes are truncated at a maximum 
magnitude |g| max = ir/2a, and the reverse Fourier trans- 
form performed with x = 0. This reverse transform can- 
not be performed analytically for arbitrary Co, C 2 and 
C4, but can be reduced by standard techniques to two 
one-dimensional integrals that can be easily performed 
numerically, 



1/ 



probe 



probe 



ft? f 1 - 3 CqC 2 + [§C C 4 - IC C 2 - Cj] s 2 + [C 2 Cj - |C C 4 + Cj + f C C 2 ] s 4 
47rai S A(s)B(s) 

/ z ext [\ Cq + (C 2 -C )u 2 
l2iraJ Q S B(s) 



(5) 



where 

A(s) = C + (3C 2 - C Q )s 2 
B(s) = l -CoC 2 + \\c C A 



2 _ 2^-0^2 



(C2C4 — C0C4 + CqC 2 ) + c 



(6) 



The itP robe equation is similar to that of ?/P robc , with /™* 
replaced by f y lrt . It is conventional to express results in 
terms of the response functions in directions parallel and 
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perpendicular to the axis of material symmetry, 

a \\ = u prob C// virt j 



U 



probe / yvirt 



(7) 



For an isotropic material, Co = 36*2 = C4 = 3G with 
G the shear modulus, and and Q can be evaluated 
analytically, 
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(8) 



matching real-space calculations [36 for an isotropic elas- 
tic body with the Poisson ratio of \ expected for affincly- 
deforming networks in 3D . 



Anisotropy and non-linearity induced by a force 
monopole 



S£ c — £ c hihjUij for an orientation n. This was then con- 
verted to a longitudinal tension t(6£) using the chosen 
single fiber model. For this calculation the extensible 
wormlikc chain model was used, which generalises the 
inextcnsible model by the inclusion of a contour modu- 
lus K [221, 



,incxt 



KIT 2 



T ' 

1+ K. 



+ ^0, 0) 



where 



t {(j>) is the dimensionless inextensible worm- 



like chain expression 



6 2ir 2 (f> 



^^^/(f)COth('^■^/(f)) — 1 



(10) 



^0 = — ^c/Mp is the natural end-to-end distance and 
k = k^Tlp the bending rigidity. The inextensible model 
is recovered in the limit K/t —¥ 00, i.e. 



If the material is deformed at any given point that 
obeys axisymmetry, and the fiber force-extension relation 
t(S£ c ) is known, then Co, G2 and C4 can be evaluated 
using ([2]), and inserted into |5]) to determine the local re- 
sponse functions a" , a 1 - ^ at that point. This procedure 
is here applied to determine both response functions at 
various distances in front of and behind an external force 
monopole f ext along the axis of symmetry, for which, 
assuming there is no spontaneous symmetry breaking in- 
duced by e.g. elastic instabilities, axisymmetry will hold. 
Even assuming hypoelasticity, it is difficult to solve the 
full problem analytically if t(S£ c ) is non-linear, hence the 
spatial response due to f ext is here determined numeri- 
cally. 

Before continuing it is important to address the as- 
sumption of affinity in this protocol. As estimated in the 
appendix, there can be significant strain gradients near 
the applied force, even when the magnitude of the force 
is somewhat low, and these can lead to a probable break- 
down of affinity at short ranges, i.e. within 1/iui of a 
lOpN force for a typical actin network. Affinity should 
not be assumed in such systems unless some compen- 
satory effect can be argued for. For other semi-flexible 
networks, affinity may be valid throughout the spatial do- 
main, in which case the procedure outlined here should be 
valid. In lieu of an alternative predictive calculation that 
does not assume affinity, we present our affinc method 
here in the understanding it may need future generalisa- 
tion to non-affine deformations. 

To perform the numerical calculations, cylindrical co- 
ordinates (r, 9, z) were employed with the z-axis aligned 
parallel to the direction of the applied force f cxt , and 
variation with 9 suppressed following the expected ax- 
ial symmetry of the solution. Displacement vectors were 
defined at regular (r, z) lattice nodes and the correspond- 
ing strains Uij were interpolated onto an interpenetrat- 
ing staggered mesh by first-order finite differencing. The 
local strains were then converted to changes in fiber 
segment length 5£ c by assuming affine deformation, i.e. 



S£ c = fu" 



(11) 



For (f> < — 1, ( 10 ) is undefined and we assume such fibers 



have buckled and no longer contribute to the mechanical 
response of the network. We employed a high contour 
stiffness, A^=500pN for vimentin-like network parame- 
ters and A'=3000pN for actin-like networks, but varying 
K made little difference to the predictions and was pri- 
marily employed to remove the singularity at u = | in 
the (inextensible) force-extension curve, aiding numerical 
convergence. 

Numerically inverting ^ and ( 10 ) gives the tension 



as a function of extension S£ c and hence strain u^j , which 
is then integrated over all orientations to give the stress 
field Gij using (fll). The equations of mechanical equi- 



librium in cylindrical coordinates with ^-variation sup- 
pressed are [31] 



= d r o~ rr + d z a rz 



= d r a rz + d z a zz + — + /fW(x) (12) 
r 

at each node, and our goal is to find the displacement 
field corresponding to global equilibrium. Here the ex- 
ternal force is localised at the origin by the function 
W(x), which integrates to one and decays to zero for 
large |x|. However, we do not use a delta function 
as this introduces sharp gradients and numerical diffi- 
culties. Instead a smooth Gaussian profile was used, 
W(x) 



V27TO- 



-(r +z )/2a ! with a = 1/im _ 



To solve (12), Newton's method was used as follows. 



The displacements at each interior node was written as a 
combined vector U, where U contains the radial and ax- 
ial components of each nodal displacement. The internal 
forces djdij were rewritten as a nodal force vector F lnt 
arranged identically to U. Matching F mt was a constant 
vector F ext giving the external force (combined with the 
weighting function) for each node. Each component of 
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F mt was determined from the relative displacements of 
nearby nodes, and was therefore a smooth function of U 
that can be expanded about a given U , giving the global 
equilibrium equation 

= F^+F^iU) = F^+F^iUoj+A^SUp+OiSU 2 ) 

(13) 

where U = Uq + <5U and the a, j3 subscripts refer to coor- 



dinates in the U and F vectors. Here, A. 



a/3 



8F^/dU a 



is a large, symmetric stiffness matrix. Neglecting the 



quadratic and higher terms in SXJ in (131 gives a linear 



equation that can be solved via matrix inversion, 



SU a 



A <*P \?l 



-ext 

p 



pint 



(Uo)] 



(14) 



The full non-linear equations are then re-linearised about 
this new point Ui = Uo + 8TJ and inverted once more, 
producing a succession of estimates U2, U3 etc. that 
obey equilibrium in their corresponding linearised sys- 
tems. Iteration continues until the largest change in any 
single component of U changes by less than some small 
threshold value, which was taken to be 10 _5 /zm here. 
Dirichlet boundary conditions corresponding to the 



known linear solution 



for a point force f were 



imposed. Varying the linear system size by a factor 
of 2 had only slight effects (roughly 1%) on the stress 
field, so our findings are not sensitive to this choice of 
boundary condition. Results presented here correspond 
to < r < 80/iTO and — 80^m < z < 80/im, with a mesh 
size of 0.4/zm in both directions. The matrix inversion 
( 14 ) was solved using the conjugate gradient method [37] , 



III. APPLICATIONS 

Here we present predictions of the calculations of 
Scc.[H]for fibers obeying the wormlike chain model. This 
model was chosen since the only alternative in common 
usage, i.e. elastic beam models [THl HOI [331 H2] , ex- 
hibit no non-linearity in their force-extension curves and 
hence no anisotropy in the hypoelastic limit. The moduli 
and response functions for a predefined prestress is con- 
sidered first, before combining these calculations with the 
spatial response field generated by a force monopole. 



aligned at an angle 8 to the axis of symmetry will be 
extended by a relative amount 



Sic ( 2n 1 • 2 

— = 7 cos 2 6 — sur 1 

lr \ 4 



(16) 



This can then be inserted into (J2J, along with a specific 
choice of t(£ c ), to determine the elastic moduli. For the 
wormlikc chain model, the tensions and extensions can be 
expressed in normalised forms by scaling the prefactor to 
"f£ p /£ c = (8£ c /£ c )(£ p /£ c ), which then fully specifies the 
magnitude of the prestrain as per (111. Fig. [2] shows the 
independent elastic coefficients for varying 7 with fixed 
v = 1/4, showing the expected stiffening along the z- 
axis with 7 > 0, with C 4 > C 2 > C . For 7 < 0, the 
order becomes Co > C2 > C4, and transverse modes now 
become stiffer due to their extension (for v > 0). Also 
note the divergence at j£ p /£ c = | which arises when the 
end-to-end distance of fiber segments aligned with the z- 
axis equals their contour length, and can extend no more 
without contour stretching. 
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FIG. 2: The independent elastic moduli Co, C2 and C4 as a 
function of the normalised prestrain magnitude yl p /£ c , for a 
Poisson ratio v = 4. Each modulus has been scaled to its 
zero prestrain values Cy fc! , which obey Cq = 3C2 = C%. The 
vertical dashed line is at y£ p /£ c = |. The thin horizontal and 
vertical solid lines are to guide the eye to the 7 = solution. 



A. Moduli for a given prestress or prestrain 

Although in many applications the perturbation will 
be an applied force, thus generating a prestress, it is also 
insightful to consider prestrains. For an axisymmetric 
prestrain, the strain tensor can be written (in Cartesian 
coordinates, with the axis of symmetry along the z-axis) 
as 



= 7diag(-z/,-2/,l) 



(15) 



where the dimensionless v parameterises the anisotropy 
(for linear response, v is simply the Poisson ratio). As- 
suming affinity, and taking v = \ [32, ,33] ), filaments 



Following Morozov and Pismen |29j , axisymmetric pre- 
stresses are parameterised by an overall magnitude To, 
and a dimensionless degree of anisotropy /?, where ft pa- 
rameterises the degree of anisotropy in a manner that 
depends on the sign of r : (3t > corresponds to a 
more positive stress in the z-direction, and /3r < to 
more positive stresses in the transverse plane. The pre- 
strain is isotropic or absent if (3tq = 0. Then the tension 
in a fiber aligned with an angle 9 to the axis of symmetry 

is \m 



t{9) = t cos 2 6 + 



1-/3 
1 + /3 



■ sm 



(17) 



G 



where t{0) > corresponds to tension and t(9) < to 
compression. For the wormlike chain model the prestress 
magnitude r can be normalised to Tq£ 2 /ktt 2 as per (11 1. 
Determining the elastic coefficients using ([2| now requires 
converting each tension r to an extension, and expand- 
ing the force-extension curve about this point, which 
is straightforward to perform numerically. The varia- 
tion of the elastic moduli with tq for a fixed anisotropy 
P = | is presented in Fig. [3J There is a clear stiffening 
of all moduli for networks under tension To > 0, with 
C 4 > C 2 > C , and a corresponding softening for r < 
with Co > C2 > C4. Also marked on this figure are two 
special tensions: tq£ 2 /kit 2 = —1, which is when fibers 
aligned with the z-axis become contracted beyond the 
range of the wormlike chain model, which we interpret as 
a buckling event beyond which the single-fiber response 
is identically zero. There is however no clear signature of 
buckling until t £ 2 c /ktt 2 = (l + /3)/(l-/3) = -2 for this P, 
when all fibers, including those oriented transversely to 
the z-axis, buckle. All linear elastic moduli equal zero at 
and below this point. 



The trends observed are consistent with the observations 
of the previous section, i.e. a stiffening for tq > with 
a" < a 1 - and both less than their zero prestress values, 
with the opposite trend for To < 0. The divergence of 
both a" and a 1 - at t £ 2 /ktt 2 = —2 (for /3 = |) cor- 
responds to the vanishing of the elastic coefficients in 
Fig. [3j There is no divergence for (3 = 1, since for this 
extreme anisotropy fibers oriented perpendicular to the 
z-axis never become prestressed. 

It is evident from Fig. [4] that, for all To, the degree 
of stiffening or softening for /3 = 1 is reduced compared 
to j3 — |, reflecting the lower net tension for a given tq 



in (jlTh . This effect is more clearly seen by plotting the 
variation with j3 for a fixed To > as in Fig. [5j demon- 
strating that the material is softer in both directions for 
P > 0, but more so in transverse directions, with the op- 
posite trend for j3 < 0. The vanishing a"' 1 - as /3 — >• — 1 
corresponds to a diverging pretension for transversely- 
aligned filaments as per (17). 
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FIG. 3: The independent elastic moduli for an axisymmetric 
prestress of normalised magnitude tq£ 2 c / 'ktv 2 , for an anisotropy 
parameter /3 = |. The moduli are scaled to their zero pre- 
stress values Cfj k i. The vertical dashed and dotted lines cor- 



respond to to£ c /ktt 



-1 and —2 respectively. The thin 



horizontal and vertical solid lines are to guide the eye to the 
to = solution. 
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FIG. 4: Response functions in parallel and perpendicular di- 
rections to the axis of symmetry versus normalised prestress 
to for P = I and P = 1 as denoted in the legend. The response 

functions are normalised to their zero prestress aj \' . 



C. Anisotropy due to a force monopole 



B. Local response in a predefined stress field 

The local displacement u due to a point force f can 
be quantified by the response functi ons a " — u z /f z 



and a = u x /f x as described in Sec. II B so higher a 
correspo nd to a softer material. Continuing with pre- 
stresses (17), a 
d5 



can be numerically evaluated using 
Examples for P = \ and p — 1 are shown in Fig. 4 



Localised perturbations to the network, in whatever 
form they take, that are of sufficient magnitude to gen- 
erate a non-linear elastic response can induce anisotropy 
in a spatially- varying manner. For simplicity we consider 
here just the case of a force monopole that can be gener- 
ated by e.g. perturbing a probe particle with an optical 
or magnetic trap, and focus attention on points directly 
behind and in front of the force, for which the axisymmet- 
ric theory developed here should be valid. This protocol 
has been recently applied to vimentin networks 38 . The 



methodology was outlined in Sec. II C 
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model (11), the affine, hypoelastic constitutive equation 
([!]) can be rewritten in terms of the normalised tension 
and extension as 



(18) 



where the normalised tension d> is the inverse of the nor- 



malised extension u in (11). Here, the normalised and 



dimensionless starred tensors are 

£ 2 



£r Uij 



2 "tJ i 



(19) 
(20) 



i.o Then the (Cartesian) force balance equations for a force 
monopole, diOij + / cxt <5(x) = 0, become 



FIG. 5: Response functions versus anisotropy /3 for fixed 
Toil/ kit 2 = 1. 



Examples of the displacement fields induced by a small 
(lpN) and large (300pN) applied forces are presented in 
Figs. [6] and [7] respectively, for material parameters repre- 
sentative of the vimentin networks of [35]. For the smaller 
force, the displacement field is symmetric behind and in 
front of the force, and indeed can be shown to obey the 
the linear solution except near the point of force applica- 
tion. By contrast, for the larger force there is a marked 
fore-aft asymmetry in which the displacements behind 
the probe exceed those in front. This is a consequence 
of the strain stiffening under tension and softening un- 
der compression that is inherent to the wormlike chain 
model. 

The non-linear force-extension curve of the wormlike 
chain model is responsible for the phenomenon of stress 
focussing, where the stresses directly behind the force 
are enhanced. As evident in Fig. [8] the normal stresses 
across a plane normal to the direction of the force, and 
directly behind it, are increased manyfold with respect 
to the equivalent linear solution. Conversely, stresses de- 
crease and fall below the linear solution at greater lateral 
distances. Since the net force on any closed surface en- 
closing the monopole must balance f cxt , an increase in 
stress in one region is expected to be compensated by 
a corresponding decrease elsewhere. What is not trivial 
here is the location of the increase, which is a consequence 
of the non-linear constitutive equation. 

Examples of ofi'- 1 are given in the inset of Fig.[9j which 
shows plots of the response functions 10/^m from the ap- 
plied force for network parameters chosen to be represen- 
tative of actin |18j and vimentin [23j . The softer vimentin 
network is clearly much more strongly affected for any 
given force / that the stiffer actin network, as expected. 
It is convenient to present these results in normalised 
forms that can be generalised to a broader range of net- 
works. To this end, note that, for the wormlike chain 



= d,u* 



/r 



(21) 



which suggests a partial normalisation for the external 
force. To make it dimensionless, we incorporate the dis- 
tance R from the applied force monopole to the point at 
which the response functions are measured, 



/ ext lg 

pKTT 2 R 2 



r 



(22) 



Replotting the response functions for the different net- 
work parameters (but the same R) against this quantity 
shows clear collapse as demonstrated in the main figure of 
Fig. [9] confirming the validity of (21 ). This means that it 
is only necessary to provide curves for different distances 
R from the applied force; once the response functions are 
normalised to their zero prestress values a"' -1 ", which are 
given in (J8| (using known expressions for G |23j). and the 
force scaled as just described, the response functions are 
fully specified. Varying R therefore generates a family of 
curves that can be used for fitting purposes to estimate 



10 



unknown parameters. To this end, we present in Fig 
curves for R = 5/im, 10/im and 20/im. The shapes of 
these curves depend on the non-linearity and we have 
not been able to derive a simple scaling with R that col- 
lapses them. Some empirical observations are mentioned 
in the discussion. 



IV. DISCUSSION 

The relevance of our findings to actual networks will 
depend on the likelihood that internal or external forces 
can be of sufficient magnitude to place the constituent 
fibers into their non-linear response regime, which is 
problem dependent. For microrheology experiments such 
forces can always (in principle) be reached, but they may 
also occur naturally. As an instructive example, tak- 
ing typical actin parameters to be l„ = 17/im and £ c = 



2.2/mi 



the unit normalised tension tq£ 2 /kit = 1, 
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10 //m 



4 xlO V m 



10" >m 



4 xlO" 4 Aim 



FIG. 6: The displacement field for an applied force of lpN 
in Euler (pre-strain) coordinates for £ p = 0.5/xm, i c — OApm 
and p — 16.25/im~ 2 . The direction of material displacements 
are given as white arrows, and the magnitude given by the 
contours as denoted by the calibration bar. The location and 
direction of the applied force is denoted by the large, open 
black arrow at the center. 
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FIG. 7: The same as Fig. [6] for an applied force of 300pN. 



when non-linearities will occur, corresponds to an actual 
tension tq rs 0.14pN. This is well within the range of 
forces capable of being generated in physiological condi- 
tions, such as by a single mysoin-II molecular motor [3], 
confirming the relevance of non-linear fiber response to 
acto-myosin mixtures [39] . It is expected that similar 
arguments will suggest the relevance of non-linear me- 
chanical fiber response in a broader range of networks 
and problems. 

Future work could aim to reduce the reliance of the 
various assumptions that were made to close the equa- 
tions, albeit at the likely expense of an increase in com- 




FIG. 8: Stress focussing for £ c — 0.55 fim, £ p = 0.5 /jm and 
p = 16.25 pm~ 2 . The normal component of the stress tensor 
parallel to the force and 10 pm behind it, as denoted in the 
inset, is plotted scaled to the corresponding linear solution. 
The 6 applied forces are shown in the legend in the same order 
as the centre of the response curves from top to bottom. 
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fit /{pK-K 2 R 2 } 



FIG. 9: The local response functions at a distance 10/im 
behind a force monopole / in parallel a" (closed symbols) 
and perpendicular a ± (open symbols) directions. Negative 
forces correspond to locations in front of the force. The re- 
sponse functions are scaled by their respective / = values, 
and the force is scaled by £ 2 /(pktv 2 R 2 ). The inset shows the 
same data plotted versus the unsealed force. Parameters were 



■ 17^im, 
0.5pm, 



2.2/im and p 
0.6/im and p - 



- 39/im for actin [T 
5/im -2 for vimentin 



and 



23 



plexity. Affinity and hypoelasticity will not always ap- 
ply, but whereas the framework for hyperelasticity {i.e. 
finite strains) is well established [271 14"0] , predictive cal- 
culations for the degree of non-affinity has thus far been 
determined largely numerically from microscopic mod- 
els. Given the prohibitive number of degrees of free- 
dom required to simulate 3D networks spanning the same 
length scales as presented here, an analytical framework 
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111 Figs.[9jand[l0 were demonstrated to collapse after scal- 
ing the external force as per ( |22| ), However, the factor 
R 2 was chosen purely to make the force dimensionless. 
It appears that, for the range of R studied, these curves 
can also be collapsed by scaling by R as f£ 2 /(R 5 / 3 pKtr 2 ), 
but only for f > 0; the / < regime remains distinct. 
This piecewise collapse, which may well be approximate, 
presumably stems from the solution to the non-linear 
elasticity problem, but we have been unable to discern 
any simple reason for the R exponent of ss 5/3, neither 
what other length scale may be included to make the 
final quantity dimensionless. Further work to elucidate 
this observation would aid in the fitting of theory to nu- 
merical data, as it would mean that only a single curve 
for one value of R would need to be evaluated to deter- 
mine the full range of axisymmetric non-linear response. 



FIG. 10: Normalised response curves for the vimentin param- 
eters of Fig[9] and R = 5pm, 10/im and 20pm, with closed 
symbols for a" and open symbols for a ± . 



Appendix: Estimate of strain gradients due to a 
point force 



is required, but so far none has yet been devised even 
for the linear response, although there has been recent 
progress |3T]. 

Finally, we remark on an observation that we do not 
yet have an explanation for. The induced response curves 



To estimate the strain gradient at a given distance from 
an applied force, we employ the known response due to 
a point monopole in a linear isotropic body. The rank-3 
tensor dkUij evaluated at a point r = rr relative to an ex- 
ternal point force f = /n is found by twice differentiating 
the displacement held [50] . 



dkUtj 



f 



167r/i(l — v) r 3 



Sijfii [Su - 3ffcf|] - 3hi [Sikfjri + Sjkhn + Sikhrj - ^fifjfkh) 
-(1 - 2v) [hj{8 ik - 3fjffe) + hi(Sj k - Sfjfk)} 



(A.l) 



where the prefactor includes the shear modulus p and the 
Poisson ratio v. To determine the characteristic magni- 
tudes, we consider all non-zero components at a point in 
front of the force f = h, and in the perpendicular direc- 
tion with n • f = 0. Explicit evaluation reveals that, of 
all these, the largest in magnitude is dfUff = f /2npr 3 
evaluated in front of the force. 

The assumption of affinity is expected to break down 
before forces / for which this largest strain gradient ex- 



ceeds the inverse crosslink length l~ x , or / ~ pr 3 i~ x . (It 
may of course break down much sooner than this, includ- 
ing in the linear regime |42j). Taking values representa- 
tive of actin networks, i.e. £ p — 17 pm, £ c = 2.2/mi and 
p = 39/im~ 2 (so p s» 25Pa) [T5], this estimate suggests a 
breakdown in affinity by the time / /r 3 ~ 10 pN/pm 3 , or 
for forces of lOpN at a range of lpm. The same calcula- 
tions for vimentin networks [23 suggests a lower force of 
just lpN for the same distance. 
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